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Extensive numerical studies have demonstrated that the two-dimensional single-band Hubbard 
model contains much of the key physics in cuprate high-temperature superconductors. However, 
there is no definitive proof that the Hubbard model truly possesses a superconducting ground-state 
or, if it does, how it depends on model parameters. To answer these longstanding questions, we 
study an extension of the Hubbard model including an infinite-range d-wave pair field term which 
precipitates a superconducting state in the d-wave channel. Using exact diagonalization on 16-site 
square clusters, we study the evolution of the ground-state as a function of the strength of the 
pairing term. This is achieved by monitoring the fidelity metric of the ground-state, as well as 
determining the ratio between the two largest eigenvalues of the d-wave pair/spin/charge density 
matrices. The calculations show a d-wave superconducting ground-state in doped clusters bracketed 
by a strong antifcrromagnotic state at half-filling controlled by the Coulomb repulsion U and a weak 
short-range checkerboard charge ordered state at larger hole-doping controlled by the next-nearest- 
neighbor hopping t' . We also demonstrate that negative t' plays an important role in facilitating 
d-wave superconductivity. 

PACS numbers: 71.10.Fd, 74.25.Dw, 03.65.Aa, 74.72.-h 



I. INTRODUCTION 

Originally proposed for the study of transition metals, 
the Hubbard model^ continues to be one of the simplest, 
and yet complex models, used to describe strongly corre- 
lated systems and possibly the cuprate high-temperature 
superconductors^'^. However, despite over 40 years of 
intensive numerical and analytical studies and over 20 
years devoted to the connection with superconductivity 
in cuprates, the repulsive Hubbard model in two dimen- 
sions has not been proven to possess a unique supercon- 
ducting (SC) ground-state in any parameter space apart 
from asympototically small coupling with a vanishingly 
small transition temperature'*'^. Promisingly, recent dy- 
namical cluster approximation calculations suggest the 
presence of finite temperature rf-wave superconductivity 
near half-filling^'^, yet the true ground-state remains in- 
accessible due to the fermion sign problem. 

Empirically, in terms of the underlying noninteracting 

band structure, the SC transition temperature is larger 
in materials possessing bonding Fermi-surfaces far from 
(tt, 0) with only a small portion nesting the antiferro- 
magnetic reciprocal lattice vector Q = (7r,7r)*"^^. This 
degree of Fermi surface bending is related to the strength 
of the hopping parameter along the diagonals of the unit 
cell - the next-nearest-neighbor hopping t'{< 0) - in the 
single-band model for the electronic dispersion. Pavarini 
et al. give reasons why a large t' (the absolute value) 
might favor superconductivity due to the minimization 



of the hybridization between cooper and apical oxygen 
orbitals*. However, most numerical studies of supercon- 
ductivity in the single-band Hubbard model do not show 
a strong dependence on t' - in fact most studies of the t-J 
and the Hubbard model indicate that negative t' hurts 
superconductivity^^"^^. Thus, this intriguing clue to the 
root cause of superconductivity remains unexplained. 

In addition, the "phase diagram" of materials in 

the cuprate family has displayed several tendencies to 
show charge ordering in doped systems, forming either 
static stripes as revealed in neutron scattering stud- 
ies of LaBaCuO and related compounds^®'^^, or local 
checkerboard- type order on the level of a few lattice spac- 
ings as revealed by scanning tunneling spectroscopy on 
SC BiaSraCaCuaOg+o. and BiaSraCuOe+x^""^^- More- 
over, recent angle-resolved photoemission^^'^'' on nearly 
optimally doped BiaSraCuOg+x has revealed a density- 
wave type of order emerging below temperatures T* 
above T^. This propensity for charge ordering in doped 
systems gives an indication that a density-wave type of 
instability may also lie in close proximity as a competing 
phase to the dominant SC order parameter, and it sug- 
gests that higher SC transition temperatures might be 
obtained as a result of controlling the emergence of the 
density-wave phase. After such a long period of study, it 
remains an open question whether this model has enough 
physics to capture the essence of the ground-state phases 
in the cuprates as a function of doping or other material 
parameters. 
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Our aim is to investigate the conjecture that the Hub- 
bard model could possess a d-wave SC ground-state, as 
well as other postulated phases that either compete or co- 
exist with d-wave superconductivity. For this, we employ 
numerical calculation in small clusters to study an exten- 
sion of the Hubbard model, including an infinite-range 
rf-wave pair field term of strength A, which precipitates 
superconductivity in the c?-wave channel. The method 
of adding a pair field term to a standard low-energy 
Hamiltonian was introduced earlier by Rigol, Shastry, 
and Haas (RSH) for the study of superconductivity in the 
t-J modeP*'^^. Our Hamiltonian guarantees a d-wave SC 
ground-state for sufficiently strong A. By gradually de- 
creasing A, we expect a quantum phase transition (QPT) 
or, in a finite-size cluster, a quantum critical crossover 
(QCC), provided the Hubbard model does not favor a 
d-wave SC ground-state. Otherwise, if no QPT or QCC 
is observed, we infer that the bare Hubbard model pos- 
sesses a d-wave SC ground-state. 

We establish a criterion for the existence of the QCC by 
monitoring the fidelity metric g (Ref. [30-33]) as a func- 
tion of A. We perform exact diagonalization (ED) calcu- 
lations on the 16B Bctts cluster^'* (16-site square cluster) 
for the two-dimensional(2D) single-band Hubbard model 
near half-filling for different values of the Coulomb re- 
pulsion U and the next-ncarest-neighbor hopping f! . We 
also look at the strongest competing order that appears 
in proximity to the SC phase. For this purpose, we study 
the QCC between li-wave SC and other phases as A grad- 
ually decreases by calculating the fidelity metric and the 
corresponding density matrix eigenvalues and their ra- 
tios, as discussed below in Eqs. (7)- (9). 

After we obtain QCC information for the hole-doped 
region of interest, we plot the phase diagram. As a non- 
traditional phase diagram, it tracks the critical strength 
of the d-wave SC term, separating the SC dominant phase 
and other phases, as a function of doping. The phase at 
A = is simply the phase of the Hubbard model with 
no additional pairing term. We emphasize that the d- 
wave SC state in the phase diagram is robust, since this 
state is favored by the pairing term added to the stan- 
dard Hubbard model, and the fidelity metric is robust in 
tracking the phase transition as the pairing term is con- 
tinuously turned off. We examine the other competing 
orders in a less robust manner, by calculating the dom- 
inant order of each phase via studying the charge and 
spin density matrices. Although finite-size scaling us- 
ing larger cluster calculations is currently unfeasible, the 
small cluster calculations could provide insights into the 
trends of the behavior of d-wave SC and other possible 
orders approaching the thermodynamic limit. As one of 
the most important results in our work, we will present 
this phase diagram in Sec. IV. 

The rest of the paper is organized as follows. In Sec. II 
we present our calculational methodology and justify the 
use of the fidelity metric to monitor the change of orders 
in ground-states. In Sec. Ill we present and discuss our 
results from finite-size cluster ED studies of the single- 



band Hubbard model for various hole-dopings and other 
model parameters. Finally, we offer a summary of our 
results and interpretations of the phase diagram in Sec. 
IV. 



II. CALCULATION METHODOLOGY 
A. Model Hamiltonian 

The single-band Hubbard model is described by the 
Hamiltonian 

HHubbard = — ^ tijc\crCja + ^ U Tlifriil, 

where (qct) creates (annihilates) an electron with spin 
a on lattice site i; tij are overlaps or "hopping" inte- 
grals restricted to the nearest-neighbor and next-nearest- 
neighbor, denoted as t{> 0) and t'{< 0); U is the on-site 
Coulomb repulsion; and nig. = cl^Cia is the number op- 
erator. 

To ensure a d-wave SC ground-state, we introduce a 
d-wave pair field term of the BCS type with a d-wave 
order parameter^^: 

Hd-wave = ~ ]^ X] ^i^i ' 

where A = A^^i+i - Ai^i+y, Aij = Ci^Cji + cj^Cii, A 
is the strength of this d-wave BCS attraction, and N is 
the total number of lattice sites. We will assume that A 
changes between and 0(1), in units of t. 

RSH have shown that for the t-J model plus a d-wave 
pair field term with strength A, a SC ground-state is 
present at some sufficiently large A in small clusters^*^'^^. 
Similar to this model, our Hamiltonian is assumed to pos- 
sess a d-wave SC ground-state when A is sufficiently large 
in any finite-size cluster, without loss of generality. Thus, 
we view A as a "continuous switch" to study the Hubbard 
model containing a d-wave SC ground-state: when A is 
large or "on" , the d-wave BCS term is dominant and the 
Hamiltonian exhibits a d-wave SC ground-state; as A ap- 
proaches zero or the "off' state, the Hamiltonian reduces 
to that of the bare Hubbard model, thus revealing the 
properties of its ground-state. This is completely analo- 
gous to the study of ground-state magnetism by slowly 
switching off an applied magnetic field. In summary, 
we begin with a large A, utilize the assumed d-wave SC 
ground-state as a starting point, and gradually decrease 
A. If a QCC does not occur, as indicated by the behavior 
of the fidelity metric defined in the following section, it 
illustrates that the Hubbard model ground-state has the 
same characteristic as the ground-state at large A: a d- 
wave SC state; if a QCC occurs, it demonstrates a change 
in the ground-state phase, indicating that the Hubbard 
model does not possess a d-wave SC ground-state. The 
nature of the non-SC state is then determined from a 
study of the density matrix, as elaborated below. 
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B. Fidelity metric 

The main point of this work, capturing a possible QCC 

between the d-wave SC state and other phases, can be ac- 
comphshed by calculating a quantity known as fidelity. 
Fidelity was first used widely in quantum information 
to provide the criterion for distinguishing between quan- 
tum states^^'"^'' and has now been utilized increasingly in 
general physics for investigating classical as well as quan- 
tum critical behavior in various systems^^"^'^''^®"''^. The 
work with fidelity on fermionic models and spin mod- 
els recently has received increased attention related to 
QPT and is powerful enough for nonstandard QPT like 
topological phase transitions, ^^'^^ where no direct exam- 
ination can be made of an order parameter. 

Given a Hamiltonian of the form H(X) = Hq+XHi, the 
overlap between the ground-states obtained with small 
changes in the parameter A 



F(A,5A) = ($o(A)|$o(A + <5A)) 



(3) 



defines the fidelity. A need not be a small parameter but 
can be tuned to lie near a critical value such that the 
ground-state wave functions change abruptly with small 
changes in A near this critical value. For small values 
of dX, |$o(A + SX)) can be expanded in the normalized 
eigenstates at A: (a = 0,1,2...); to second order 

in 5X the term proportional to the unperturbed ground- 
state is 



($o(A)|$o(A + 5A)) 



2 ^ 



\ {<p^{xm\<fo{m' 

[Eo{X) - E^{X)]^ 

(4) 

The fidelity has been shown to be an extensive quantity 
that scales with the system size in other Hamiltonians 
similar to our model^^'"^^. So, following RSH, we intro- 
duce the fidelity metric^°"^^ 



g{X,SX) = i2/N){l- F{X,SX))/6X^ 



(5) 



(equivalent to the fidelity susceptibility defined in other 
contexts"^^'**^). , whose limit as (5A — >■ is well defined 
away from the critical points when perturbation theory 
applies 



lim g(X, SX) 



-y 



N 



|(<l>»(A)|gi|$o(A))|^ 
[Eo{X) - E„{X)]^ 



(6) 



From Eqs. (5) and (6), the fidelity metric is in principle 
an intensive quantity that does not scale with the number 
of sites and measures how significantly the ground-state 
wave function changes as the parameter A changes. At 
the critical point where perturbation theory fails, the nor- 
malized ground-state wave functions differ significantly 
with only a small change in the parameter A, which re- 
sults in a significant decrease of the overlap of these two 
wave functions or the fidelity at that point. Thus, we 
expect the fidelity metric to evolve smoothly away from 
critical points but to undergo a dramatic increase (or 



drop in other cases'^") at them. By approaching a quan- 
tum critical point from either direction, the fidelity met- 
ric can be a powerful tool to signal a QPT/QCC. 

In our numerical calculations, we find that the calcu- 
lated fidelity metric is independent of the value of SX 
when SX is between 0(10^'') and O(10~^), which is con- 
sistent with its intensive behavior; however, the fidelity 
metric could behave substantially diff'erently when SX is 
smaller and larger, due to numerical errors and lack of 
accuracy, respectively. Thus, in what follows, we use 
SX = 3 X 10~^ to balance these two concerns. 



C. d-wave pair density matrix 

While the fidelity metric captures the existence of a 
QPT/QCC, we would also like to investigate quantities 
directly related to the order parameter of the ground- 
states. The d-wave pair density matrix^^'^^'^* can be 
defined as 



pd- 
ij 



{^o\dIDj\^o) 



(7) 



One important quantity is the ratio between the 

largest and next largest eigenvalues of the d- 
wave pair density matrix, denoted as Rd-wave = 
^^(^pd-wavey^^(^pd-n,ave^_ The study of the ratio 

Rd-wave was first suggested by RSH^*^'^^. To understand 
its significance, we note that based on the work of Pen- 
rose and Onsager^^, particle condensation into some well- 
defined ordered state leads to characteristic behavior in 
the corresponding density matrix: the largest eigenvalue 
scales to a finite value greater than 0(1) while the other 
eigenvalues are ~ 0{1). Further, Yang^^ showed that the 
largest eigenvalue actually scales linearly with the total 
number of particles. Thus, Rd-wave ~ "^^N + where 
^ is the order parameter and $ is sub-linear in the sys- 
tem size. In other words, the ratio Rd-wave, having an 
intrinsic correlation to this rf-wave SC order parameter, 
can be used directly to investigate d-wave SC order. A 
similar analysis follows for other (non-SC) orderings. 



D. Other quantities of interest 

Besides superconductivity, both charge and spin orders 
have been reported in cuprate high-temperature super- 
conductors, with the undoped parent compounds known 
to be antiferromagnetic Mott insulators. Moreover, the 
2D Hubbard model for sufficiently strong Coulomb repul- 
sion U possesses a well-known strong antiferromagnetic 
order at half-filling. Thus, we study charge/spin order- 
ing or the tendency toward ordering by calculating the 
charge and spin density matrices: 



(8) 



(9) 
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where a^^ = &^Ci-^ - 0:^04 and pi = c^^ai + cJ^Q;. 

We investigate both i?,p,„ = Ai{P^'p")/A2{P-f") and 

Rcharge = Ai (ff^''°'''''')/A2 (P^''"''^'') , as we do for the d- 
wave pair density matrix. We also examine the specific 
matrix elements of the density matrix associated with the 
correlation functions. By calculating the above quan- 
tities, we investigate the strongest order that exists in 
proximity to the d-wave SC order in the small cluster, to 
provide insights into the behavior of competing or coex- 
isting states in the thermodynamic limit. 



III. RESULTS 



Before describing our numerics, we clarify our view- 
point. Based on the experience in RSH, we look at the 
fidelity metric g, together with the ratio R of the largest 
to next largest eigenvalues of the density matrix with 
suitable symmetry. It is a combined look at these vari- 
ables g and R as functions of A that enables us to draw 
conclusions about the nature of ordering. We note that a 
peak in g usually signifies a crossover between two phases, 
and a jump in g signifies a level crossing. We denote by 
A* the location of the peak in g{X), and it is seen below 
to play an important role in the analysis. The ratio R 
has a magnitude ~ 1, when the symmetry of the order 
parameter is incompatible with the ground-state. It is 
large (3-10 in the t-J model studied in RSH) when the 
symmetry is compatible, for small systems. In the limit 
of large systems, which is currently unattainable, it is 
not possible to assert that a quantum critical point ex- 
ists, but looking at g and R together gives us a reasonable 
picture of the evolution of the ground-state. 

We perform the calculations with the ED technique in 
a finite-size cluster using the Parallel ARnoldi PACKage 
(PARPACK)""^ based on the Implicitly Restarted Arnoldi 
method. This approach is efficient in obtaining a set of 
properly orthonormalized eigenvectors and is more nu- 
merically stable compared with standard Lanczos meth- 
ods. In our ED problem, the size of the Hilbert space 
grows exponentially with the increase of the cluster size. 
Due to computational time considerations, we perform 
our calculations for the 16B Betts cluster'^'' with periodic 
boundary conditions. Moreover, the 16B Betts cluster 
provides access to momentum (7r/2,7r/2), which is im- 
portant in catching low energy excitations at half-filling 
and low dopings. Considering symmetries, including the 
invariance of the total number of electrons and total mo- 
mentum, and fixing ourselves in the space with total 
spin Sz = 0, the Hilbert space size of the Hamiltonian 
is diminished to ~ 10^ so that the ground-state eigen- 
vector and eigenvalue could be calculated by ED using 
PARPACK in a reasonable time. 
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FIG. 1; The fidelity metric g as a function of d-wave pair field 
strength A (in units of t) at half-filling. Coulomb U {6t, 8t, lOt) 
and t'(0, — 0.15i, — Q.3i) are indicated. The small magnitude of 
g indicates the incompatibility of d-wave order at half-filling. 



A. Antiferromagnetic Order at Half-filling 

The fidelity metric g at half-filling for various values 
of U and t' (measured in units of t) is shown in Fig. 1. 
g approaches a small value as A — > 0, and 5 — as 
A increases beyond the peak value. These two regions 
correspond to two different phases: the rf-wave SC state 
at large A and the competing state at small A. A* is 
larger than 1, implying that this competing state exhibits 
strong order which is even stronger as U increases, since 
it requires a large d-wave A to conquer the competing 
phase. 

To obtain a full picture, we also calculate both Rd-wave 
and Rspin as functions of A with the results shown in 
Fig. 2. First, Rd-wave reaches 1 as A — )■ 0, which confirms 
the incompatibility of d-wave SC order in the Hubbard 
model at half-filling. Second, the large value of Rspin at 
A — > is consistent with the well known fact that the 2D 
Hubbard model at half-filling exhibits strong antiferro- 
magnetic order. The values of Rspin at A — >■ increase 
for increasing U, demonstrating stronger antiferromag- 
netic order for larger values of U. Further examining the 
spin correlation functions (not shown) reveals a (7r,7r) 
antiferromagnetic order. Third, with the increase in A, 
Rspin decreases and Rd-wave increases, with a crossover 
at a certain A. The crossover of Rspin and Rd-wave in 
Fig. 2 moves to higher A at larger U, consistent with the 
A* from the fidelity metric g in Fig. 1. Based on these 
points, the fidelity g in Fig. 1 is now more clearly un- 
derstandable. The small peak at A* between 1 to 2.5 for 
various parameters arises from forcing d-wave SC order 
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FIG. 2: The ratio of the largest to second largest eigenvalue 
of the rf-wave pair density matrix Rd-mave and that of the 
spin density matrix Rspin as functions of d-wave pair field 
strength A (in units of t) at half-filling. Coulomb U (6i , 8t, lOt) 
and f'(0, — O.lSt, — 0.3f) are indicated. 



by the applied SC term in Eq. (2). The peak should grow 
with increasing cluster size to indicate a true QPT, and 
thus for our cluster the peak implies a crossover (QCC) 
rather than a true QPT. In a strict sense we are thus 
studying QCC rather than quantum critical points. 

We next discuss the effects of t' . Figure 1 shows that 
the QCC moves to lower values of A as the magnitude of 
t' increases. Moreover, Fig. 2 demonstrates that negative 
t' enhances only slightly the d-wave SC order measured 
by Rd-wave while substantially suppressing the antiferro- 
magnetic order as reflected in Rspin- This indicates that 
while negative t' may hurt the antiferromagnetic order, 
it does not necessarily lead to an enhancement of super- 

and 



conductivity at half filling. These effects in Rgpm 
Rd-wave do uot depend critically on the value of V 
this range of parameters. 
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FIG. 3: The fidelity metric <; as a function of d-wave pair 
field strength A (in units oft) at 12.5% hole-doping. Coulomb 
U{6t, 8t, lot) and t'{0, -O.I5t, -O.St) are indicated. 



in the momentum sector of the ground-state. This type 
of sudden change for these sets of parameters only occurs 
for this specific cluster and should not occur in the ther- 
modynamic limit. The trends in Rd-wave and Rspin as 
functions of A at various values of U and t' are also simi- 
lar to the results at half-filling shown in Fig. 2 — d-wave 
SC and spin order depend mainly on U and slightly on t' . 
As seen in Fig. 4, Rspin is larger than 1 as A 0; how- 
ever, it is smaller here than at half-filling, implying that 
the competing state is also antiferromagnetic but weaker 
than that at half-filling. The sudden jump in Rspin for 
t' — — 0.3t in Fig. 4 occurred for the same reason as that 
in g, as we have discussed. 

In contrast to half-filling, at this doping negative t' 
suppresses antiferromagnetism and slightly enhances su- 
perconductivity, as shown in the larger value of Rd-wave 
and the smaller value of Rspin for negative t' at A '-^ — 2 
in Fig. (4). These trends exist for all the values of U in 
our calculations. 



B. Hole-doping near Half- filling 

Figure 3 shows results at 12.5% hole-doping that fol- 
low a similar trend as half-filling: as A — 0, g approaches 
a small finite value; as A increases, g increases to a max- 
imum at some A* and then decreases monotonically to- 
ward zero with additional increase in A. This behavior 
indicates that the ground-state of the Hubbard model is 
not a d-wave SC for this cluster size or parameter regime. 
The dependence of A* on U and t' is also similar to the re- 
sults at half-filling, except for the discontinuities in g for 
t' = — 0.3t at A ~ 1.0, stemming from an abrupt change 



C. SC state at 25% hole-doping 

Figure 5 shows the fidelity metric calculations at 25% 
hole-doping, with a log-linear scale employed to catch 
clearly the transition as A — ?► 0. As seen in Fig. 5, A* 
moves close to zero or even negative values (see also 
Fig. 13), suggesting that the state remains in the same 
phase as A decreases toward zero. In Fig. 6, the robust 
value of Rd-wave^ 10 for A ^ 3 is an indicator of d-wave 
SC order; Rd-wave persists to ~ 2 (but is not as small as 
~ 1) for A ~ 0. This strongly suggests the existence of 
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FIG. 4: The ratios of the largest to second largest eigenvalue 
of the d-wave pair density matrix Rd-wave and that of the spin 
density matrix Rapin as functions of d-wave pair field strength 
A (in units of t) at 12.5% hole-doping. Coulomb U = 8t and 
t'(0, -O.lSt, -0.3t) are indicated. 



d-wave SC order in the 2D single-band Hubbard model 
even without the d-wave pair field. We note that this is 
supported by recent calculations based on infinite pro- 
jected entangled-pair states in the t-J model which yield 
the strongest SC pairing around 25% doping.'''' 

There is growing evidence that t' plays an impor- 
tant role in the critical temperature and other proper- 
ties of SC materials^"-' Our calculations support 
the statement that negative t' favors a d-wave SC state, 
indicated by the fidelity metric calculations shown in 
Fig. 5, in which negative <' gives a rapid increase to 
the fidelity metric as A goes to zero and yields A* — > 
with increasing magnitude of f , and also demonstrated 
by the d-wave pair density matrix calculations, in which 
Rd~wave (at A = 0) has a larger value as the magnitude 
of t' increases. Notice there are sudden changes of g at 
t' = OandU = 6t, 8t or lOt, only near A — ?> 0, which do 
not originate from a change of the momentum sector of 
the ground-state as encountered in previous cases. These 
sudden changes are seen in Fig. 6 for Rd-wave as well, 
where Rd-wave jumps to 1 as A —> 0, where the eigen- 
states of the d-wave pair density matrix with the largest 
eigenvalue are two-fold degenerate. 



D. Hole-doping far from Half- filling 

In the highly doped region (37.5 and 50% hole-doping), 
fidelity metric calculations display positive A*, indicat- 
ing a ground-state incompatibility with d-wave pairing 
as A — )■ 0, and give a A* that depends strongly on t': 
at 37.5% hole-doping [Fig. 7], A*- 0.25,0.5 and 0.75 at 
t' = -0.3i, -0.15i and 0; at 50% hole-doping [Fig. 8], 
A*- 0.25,0.4 and 0.6 at t' = -0.3t, -0.15i and 0. Since 
A* here is one order of magnitude smaller than the A* 
in the weakly doped region, the competing state here 
is more easily suppressed by the d-wave pair field when 
compared with the strong antiferromagnetic order exhib- 
ited near half-filling. 

Calculations of the charge density matrix show 
that Rcharge IS always greater than 1 and does 
not change much as a function of A, so tracking 




FIG. 5: The fidelity metric gr as a function of d-wave pair 
field strength A (in units of t) at 25% hole-doping. Coulomb 
U{6t,8t,Wt) and t'(0, -0.15t, -0.3t) are indicated. A semi- 
log scale is employed to highlight the behaviors around A = 0. 
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FIG. 6: The ratio of the largest to second largest eigenvalue 
of the d-wave pair density matrix Rd-wave as a function of d- 
wave pair field strength A (in units of t) at 25% hole-doping. 
Coulomb U{6t,8t,lQt) and t'(0, -0.15t, -0.3f) are indicated. 
The behaviors near A = are highlighted in the insets. 



Rcharge by itsclf might be insufficient to determine 
the change of charge order. We then focus on the 
charge correlation functions (CCFs) — the elements in 



the charge density matrix P'^ 



^charge 



after being normalized. We investigate the 

normalized nearest-neighbor CCF (NNCCF) 
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FIG. 7: The fidelity metric <; as a function of d-wave pair 
field strength A (in units of t) at 37.5% hole-doping. Coulomb 
f/(6t, 8t, lot) and t'{Q, -0.15t, -0.3t) are indicated. 
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1 

FIG. 8: The fidelity metric <; as a function of d-wave pair 
field strength A (in units of t) at 50% hole-doping. Coulomb 
U{Gt, St, lot) and t'(0, -0.15t, -0.3t) are indicated. 



($n|/OrPr+£|$o)/($o|PrPr|$o), the second-NNCCF 
2NNCCF {^Q\prPr+x+y\^Q)/{^o\PrPr\^o) and the third- 
NNCCF 3NNCCF ($o|PrPr+2i|$o)/(*o|PrPr|$o), and 
we observe how they evolve as functions of A. Figures 9 
through 11 show the normalized CCFs and Rd-wave at 
37.5% doping for U = lOt, 50% doping for U = IQt and 
50% doping for U = 6t. We observe a QCC between 



0.60 
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FIG. 9: The CCFs (NNCCF, 2NNCCF and 3NNCCF) and 
the ratio of the largest to second largest eigenvalue of the d- 
wave pair density matrix Rd-wave as functions of d-wave pair 
field strength A (in units oft) at 37.5% hole-doping. Coulomb 
U = lot and t'(0, -0.15t, -0.3t) are indicated. 




FIG. 10: The CCFs (NNCCF, 2NNCCF and 3NNCCF) and 
the ratio of the largest to second largest eigenvalue of the d- 
wave pair density matrix Rd-wave as functions of d-wave pair 
field strength A (in units of t) at 50% hole-doping. Coulomb 
U = lot and t'(0, -0.15t, -0.3t) are indicated. 
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FIG. 11: The CCFs (NNCCF, 2NNCCF and 3NNCCF) and 
the ratio of the largest to second largest eigenvalue of the d- 
wave pair density matrix Rd-wave as functions of d-wave pair 
field strength A (in units of t) at 50% hole-doping. Coulomb 
U = 6t and t'{0, ~0.15t, -0.3t) are indicated. 
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FIG. 12; The density of states for 16-site single-band Hub- 
bard model without d-wave pair field at 37.5% hole-doping 
a,t U = 8t and t' {0, —0.15t, —0.3t), energy is in units of t. 
Energy zero represents the ground-state energy. Dash line 
and solid line represent the occupied and unoccupied states 
respectively. 



the d-wave SC state and some charge ordered state by 
the trends in Rd-wave and the CCFs as A decreases 
toward {Rd-wave and NNCCF drop, and the other 
CCFs increase). This is a QCC from the d-wave SC 
state to a checkerboard charge ordered state. Careful 
examinations of CCFs at A = show that this charge 
order for the bare Hubbard model is short-ranged or 
local, since CCFs at longer distances (up to the half of 
the diagonal of this cluster - the largest possible distance 
for CCFs) do not exhibit distinctive charge order (not 
shown). Therefore, our calculation shows a short-range 
checkerboard charge ordered state at large hole-doping 
in the 2D Hubbard model. The evolution from the local 
checkerboard charge order to the d-wave SC order, as 
A increases, originates from the fact that Cooper pair 
formation favors the nearest-neighbor pairing or charge 
correlation, which destroys the checkerboard charge 
order where the nearest-neighbor charge correlation is 
small. The transition at 50% doping (see Fig. 10) is 
similar to that at 37.5% doping. One interesting finding 
is that this QCC is dependent not on Coulomb U but 
rather on i', demonstrated by the comparison between 
the calculations for different values of f/, as shown in 
Fig. 10 and Fig. 11. At this doping, electron density is 
low and double occupancy is less relevant, so Coulomb 
U has little effect on the properties of the ground-state. 

Our results show that as we turn on A and observe 
how the local charge ordered state competes with d-wave 
superconductivity, negative t' seems to only slightly sup- 



press the charge order but rather largely reinforce the su- 
perconductivity. This is indicated by comparing Rd-wave 
at different t' , where Rd-wave at a large magnitude of f 
grows much faster than at a small magnitude of t' as A 
is turned on. 

The way negative t' enforces the d-wave superconduc- 
tivity might be related to the suppression of antiferro- 
magnetic order near half-filling, as discussed in Sec. III. A 
and III.B. We also note that the gathering of density of 
states in proximity of the Fermi level at negative t' for 
our parameter range may also help form SC order. To 
justify this point, we calculate the density of states for 
the pure Hubbard model at 37.5% hole-doping, as shown 
in Fig. 12. This is done by summing over the spectral 
functions at all momentum points from our 16-site cal- 
culations. The dashed lines show the removal spectra 
or the occupied states, and the solid lines show the addi- 
tion spectra or the unoccupied states. A more negative f 
facilitates a gathering of density of states in both the re- 
moval and addition spectra at the proximity of the Fermi 
level (which can be defined as the crossover of the dashed 
line and the solid line) . The gathering of density of states 
near the Fermi level helps formation of superconductivity 
by providing more weight to mix particles and holes to 
form the condensate. Particularly at high hole-doping, 
this trend prevails over any effect that i' may have on 
antiferromagnetic ordering or charge density wave ten- 
dencies. As the density of states at the Fermi level is 
proportional to the SC susceptibility, our calculations in- 
dicate a stronger susceptibility for negative t' in the high 
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hole-doping regime for our parameter range. 



IV. PHASE DIAGRAM, SUMMARY AND 
CONCLUSION 

Based on the above calculations at various dopings, 
we obtain the Hubbard model phase diagram tracking 
A* as a function of doping, as shown in Fig. 13. Com- 
bined with the d-wave pair, charge, and spin density ma- 
trix calculations, several conclusions are obtained: (i) At 
half-filling and low hole-doping, the Hubbard model ex- 
hibits a strong antiferromagnetic order with A* depend- 
ing primarily on Coulomb U; (ii) for high hole-doping, 
we observe short-range checkerboard charge order with 
A* exclusively determined by t'; (iii) in our cluster at 
25% hole-doping. A* moves close to zero or even negative 
values, demonstrating a d-wave SC state in the pure Hub- 
bard model in this cluster, supported by other numerical 
studies''^; (iv) negative t' favors d-wave SC order, where 
A* is smaller at larger absolute value of t' and where 
Rd~wave has a larger value at A = for larger absolute 
value of t'; and (v) the antiferromagnetism at half-filling 
and low doping is much stronger than the short-range 
checkerboard charge order for high hole-dopings, revealed 
by a much larger A* in the former case. 

The phase diagram has several compelling features 
that indicate that the Hubbard model contains rich 
physics that may underlie an understanding of different 
competing phases in the cuprates. The various phases 
that our calculations suggest follow trends similar to 
those from recent density matrix renormalization group 
calculations on the infinite U Hubbard model^" and in- 
finite projected entangled-pair states calculations in the 
t-J model**^. Clearly, our results show the prevalence 
of antiferromagnetic order around half-filling and that 
negative t' only slightly assists superconductivity while 
largely suppressing antiferromagnetism. Here, the dom- 
inant energy scale is the Hubbard U. At large dop- 
ing, a short-range checkerboard charge order emerges. 
This short-range checkerboard charge order is more eas- 
ily destabilized by the d-wave pair field than antiferro- 
magnetism, giving smaller values of A* at large hole- 
dopings. Here, negative t' largely facilitates supercon- 
ductivity and slightly suppresses checkerboard charge or- 
der. In between, it appears that superconductivity may 
emerge as a consequence of the battle between antiferro- 
magnetism controlled by U and density of states and pos- 
sibly checkerboard charge order controlled by t' . There- 
fore, it is tempting to intuit that the physics of t' may 
not be reflective of orbital content (such as apical oxy- 
gen content) or a destabilization of the Zhang- Rice sin- 
glet but more connected to a destabilization of a com- 
peting spin/charge order and density of states effects at 
low and high hole-doping and a reinforcement of the SC 
order at low and high hole-doping. Thus, it would be of 
great interest to explore larger cluster sizes to examine 
the changes of the phase diagram as a function of cluster 
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FIG. 13: Phase diagram of A* (in units of t) vs. hole-doping 
rate at (7 = 6t,8t,Wt and t' = 0, -0.15t, -0.3t. The 2D 
single-band Hubbard model favors a d-wave SC phase at 25% 
hole-doping, an antiferromagnetic phase at low doping and 
local checkerboard charge order at high doping. Strong an- 
tiferromagnetic order near half-fiUing is mainly dependent 
on Coulomb U, while the strength of the local checkerboard 
charge order far from half-filling is exclusively determined by 
t'. 

size. 

In summary, ED techniques have been employed to 
study the d-wave SC phase of the 2D single-band Hub- 
bard model. The fidelity metric g, the d-wave pair and 
charge/spin density matrices, including Rd-wave, Rspin, 
and CCFs, have been calculated for the single-band Hub- 
bard model plus a d-wave pair field term at various hole- 
dopings and parameters, including Coulomb ?7, the next- 
nearest-neighbor hopping i', and the d-wave SC pair field 
strength A. Calculations show that at 25% hole-doping 
the 2D single-band Hubbard model possesses a d-wave SC 
ground-state. Based on these results, we further summa- 
rize the behavior of A* as a function of doping in the 
phase diagram. Moreover, the phase diagram shows that 
d-wave SC order is enhanced by negative i', which is 
consistent with experiments®^"'^ ^. Finally, to obtain the 
d-wave SC order parameter, finite-size scaling analysis 
could be done, provided that ED calculations on larger 
systems become feasible in the future. 
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